The QCD phase diagram according to the center group 
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We study an effective theory for QCD at finite temperature and density which contains the 
leading center symmetric and center symmetry breaking terms. The effective theory is studied in a 
flux representation where the complex phase problem is absent and the model becomes accessible to 
Monte Carlo techniques also at finite chemical potential. We simulate the system using a generalized 
Prokof 'ev-Svistunov worm algorithm and compare the results to a low temperature expansion. The 
phase diagram is determined as a function of temperature, chemical potential and quark mass. 
Shape and quark mass dependence of the phase boundaries are as expected for QCD. The transition 
into the deconfined phase is smooth throughout, without any discontinuities or critical points. 
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Obtaining a deeper understanding of the QCD phase di- 
agram wiU be one of the central goals of particle physics 
in the coming years. With running and upcoming ex- 
periments that drive this development, also the theoret- 
ical side is challenged to improve our understanding of 
the QCD phase structure. Analyzing phase transitions 
is clearly a non-perturbative problem and suitable tech- 
niques have to be applied. For vanishing chemical po- 
tential, lattice QCD is a powerful method that provides 
reliable quantitative information on the QCD finite tem- 
perature transition. However, for non-vanishing density 
the notorious complex phase problem so far limits nu- 
merical lattice QCD studies to painfully small volumes. 

For quenched QCD, where the quark contributions to 
the path integral are neglected, the deconfinement tran- 
sition is related to the center group Z3 of SU(3), which 
is a symmetry in the confined low temperature phase, 
while it is broken spontaneously above the deconfinement 
temperature [ij . When one couples the dynamics of the 
quark fields, the center symmetry is broken explicitly by 
the fermion determinant. This explicit breaking overlays 
the spontaneous breaking of the quenched theory. How- 
ever, as for spin systems, one may expect that also for 
QCD the underlying symmetry still governs parts of the 
dynamics of the full theory, e.g., via the center properties 
of canonical determinants 0. 

In order to study the role of center symmetry for the 
QCD phase diagram, we here analyze an effective theory 
which contains the leading center symmetric and center 
symmetry breaking terms. This effective center theory 
can be mapped exactly to a flux representation Q , where 
the complex phase problem is absent. So far the model 
was studied only in a very limited parameter range ■ 
Here we apply a generalization of the worm algorithm 
which allows us to efficiently explore the full range of 
temperatures and chemical potential values. This consti- 
tutes one of the very few examples where a QCD related 
complex phase problem can be solved. We study the 
phase diagram of the effective center theory and analyze 
what role the center degrees of freedom of QCD play for 
the phase structure of hot and dense matter. 



Eflfective center theory and flux representation 

The effective center theory is defined by the action Q 



s[p] = -E [p^p-^ 
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The dynamical degrees of freedom are the center vari- 
ables e = {1,6*277/3^ g-j2,r/3} ^J^g gj^gg ^ of 

a 3-dimensional cubic lattice. The partition function 
is a sum over all configurations of the center variables, 
Z = Y.{p} exp(-5'[P]). The first term of the action ^ 
is a nearest neighbor interaction which is invariant un- 
der global center transformations (P^ — >■ zP^), where all 
variables are transformed with a center element z e Z3. 
The form of this term may be obtained from a strong cou- 
pling expansion of the effective action for the Polaykov 
loop which in quenched QCD is the order parameter for 
center symmetry and for confinement. The variables Px 
take over the role of the local Polyakov loops. Although 
in full QCD center symmetry is broken explicitly by the 
quarks, the Polyakov loop is still used to monitor confine- 
ment properties and to determine the crossover temper- 
ature (see, e.g., [1, [l^l)- The strong coupling expansion 
also identifies the parameter r as an increasing function 
of the temperature T of the underlying lattice QCD the- 
ory, and for brevity we refer to r as temperature. 

The second term of the effective action may be ob- 
tained from a hopping (i.e., large quark mass) expansion 
of the fermion determinant and contains the leading cen- 
ter symmetry breaking contributions. The parameters 
77 = Ke^ , rj ^ Ke~^ are related to the chemical poten- 
tial /X. The hopping expansion shows that k, ~ Nj h{m), 
where Nf is the number of flavors and h{m) is a function 
of the QCD quark mass m which decreases with increas- 
ing TO. We refer to n as the inverse mass parameter. 

For vanishing external field, k = 0, the model re- 
duces to the 3-state Potts model, which is known to 
have a first order transition at r = 0.183522(3) 
For small external field k and vanishing fi the first or- 
der transition persists, giving rise to a short first order 
line which terminates at a critical endpoint at (r, k) = 
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(0.183127(7), 0.00026(3)) For small non-zero ^l the 
system has been analyzed with techniques based on the 
Swendsen-Wang cluster algorithm [1], with reweighting 
Q and with imaginary ^ [6|. Within the flux representa- 
tion ^3] local Metropolis updates were also used [3, .4, X]- 
It has been demonstrated that turning on the chemical 
potential softens the transition and shifts the critical end- 
point towards smaller values of k. So far no simulations 
were done in the parameter region where the complex 
hase problem of the formulation ((IJ becomes severe (see 
for a discussion of that regime). 

The flux representation solves the complex phase 
problem. We briefly summarize it to discuss our observ- 
ables and conventions: The Boltzmann factors for the 
nearest neighbor terms of ([T]) can be rewritten as 



.^[■Pxp:+^+c.c.] 



(2) 
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The sum on the right hand side is over flux variables 
bx,u G {—1,0,-1-1} attached to the links of the lattice. 
The constants C and B depend on the temperature r via 
C = (e2^ + 2e-^)/3 and B ^ (e^^ - e~^)/3C. Similarly 
one expresses the center symmetry breaking terms as 



(3) 



Sx = -1 



where we sum over monomer variables s^ G {—1,0, +1} 
attached to the sites x. It is straightforward to work out 
the monomer weights Mg for s = —1,0,-1-1, 



(4) 



The weights Mg turn out to be non-negative. 

Inserting ^ and ^ into the partition sum gives rise 
to a complete factorization of the dependence on the dy- 
namical variables Px and the sum over all configurations 
can be performed. One ends up with a new form for the 
partition sum (we drop an irrelevant overall constant). 



Z=Y^W{b,s) \{T{y.[bx,u - bx-pA + 



(5) 



{b,s} 



The partition function is now a sum over configurations 
{6, s} of the flux and monomer variables. Each con- 
figuration comes with a real non-negative weight factor 
W(b,s) = {]\x.uB^''''-"^)^Y\x ^sj- For every link with 
non-vanishing flux, i.e., bx,u = il, a factor B is taken 
into account. Sites x contribute with factors Mg^ ac- 
cording to their monomer values Sx G { — 1,0,+1}. 

The flux/monomer configurations are subject to con- 
straints: In ([5]) T{n) is the triality function defined as 
T{n) — (5„jjjQ(j3 Q. The constraints enforce that at ev- 
ery site X the total fiux from both, flux variables bx^u and 
monomers s^;, is a multiple of 3. In the flux form ([5]) the 
partition sum contains only real and non-negative con- 
tributions and thus the complex phase problem is solved. 



In this letter we focus on bulk observables such as the 
order parameter P and the corresponding susceptibility, 
which both are obtained as derivatives of the free energy, 
(P) ~ dhiZ/drj and xp = \nZ/drf' . In a similar way 
one obtains the internal energy U and the heat capacity 
C. For an efficient evaluation the identities dM+i/dr] — 
Mo, dMo/drj = M_.i and dM_i/dr] = M+i are useful. 
In the end all our observables are expressed in terms of 
the total flux and monomer numbers and their moments. 

Simulation with the worm algorithm 

Having established the flux representation where the 
complex phase problem is solved, we now must flnd 
a suitable algorithm for an efficient Monte Carlo up- 
date. We here use a generalized form of the Prokof'ev- 
Svistunov worm algorithm Q : The worm starts at a ran- 
domly chosen site and moves along links until it returns 
to the starting point where it terminates. We allow for 
two different moves of our worm: A) The worm ran- 
domly chooses a new direction at a site and changes the 
fiux at the corresponding link by ±1. B) The worm de- 
cides to change a monomer variable by ± 1 and then ran- 
domly hops to another site where the monomer variable 
is changed by =pl. The moves are offered with equal prob- 
ability, produce only configurations that are compatible 
with the constraint, and lead to an ergodic algorithm. 
The Metropolis acceptance probabilities are pA = B^'' 
when changing a fiux variable b by an amount of Sb (Move 
A) , and pb = Mg' /Mg for changing a monomer variable 
from s to s' (Move B). A more complete account of the 
algorithm and its implementation will be given elsewhere. 

We generated ensembles for lattice sizes 36'^, 48'^, 64'^ 
and 72'^. For the inverse mass parameter we used k — 
0.1, K = 0.01, K = 0.005 and k = 0.001. The evaluation of 
our observables (P), xpi U and C is based on up to 10^ 
configurations, separated by 10 worms for decorrelation. 
Autocorrelation times were determined and used in the 
estimate for the statistical errors. Finite volume effects 
were analzed by comparing the different system sizes and 
are negligible for our final results. 

The results from the new worm algorithm were checked 
using several strategies. For vanishing inverse mass pa- 
rameter the known results [llj for the 3-state Potts model 
with external magnetic field were reproduced. For small 
values of r and arbitrary k and /i we used low tempera- 
ture expansion techniques to determine the power series 
for the partition sum, taking into account all terms up 
to . For small r we found excellent agreement between 
the Monte Carlo results and the perturbative series (see 
Fig. [5]) . Finally, for all our production and analysis codes 
two independent programs were written for cross checks. 

Results from the Monte Carlo calculation 

We begin the discussion of our results with the order 
parameter (P) , which - as discussed above - is identified 
with the Polyakov loop of QCD. In Fig. [1] we show the 
results for {P)/V as a function of r and n for our 36'^ 
ensembles at k = 0.01. For roughly 450 points in the 
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FIG. 1: The order parameter (P) /V as a function of tem- 
perature r and chemical potential /i from our 36^, K = 0.01 
ensembles. Near the rear left corner no data were computed. 
(P) /V has values close to 1 there. For small r and p the order 
parameter approaches 0. 

T-/i plane the values of {P)/V were evaluated and used 
for the 3-D plot. In the rear left corner and for small r 
and /i no data were computed. (P) /V is expected to be 
close to 1 in the rear left corner. For small r and /x there 
is a sizable region where the expectation value {P)/V is 
small and center symmetry is broken only very mildly. 
Transferring this finding from the effective center theory 
to QCD implies that for small temperature and density, 
matter is confined. When r or /i are increased, the system 
undergoes a change and (P) /V reaches values close to 1. 
For QCD this implies that both temperature and fi may 
be used to drive the system into the deconfined phase 
characterized by a large Polyakov loop. 

The next step in the analysis of the phase diagram is to 
identify the phase boundary. For that purpose we studied 
the susceptibility xp and the heat capacity C as a func- 
tion of /X at fixed r (symbols with horizontal error bars 
in Figs. 2 and 3) or as a function of r at fixed (vertical 
error bars) and determined the position of the maximum: 
We fitted the data for Xp and C near the maxima with a 
parabola and obtained the position of the maximum as 
one of the fit parameters. The corresponding statistical 
error was computed with the jackknife method. In Fig. [5] 
we show the positions of the maxima of xp in the r-/i 
plane. We compare the results for 4 values of the inverse 
mass parameter k and connect the data at the same k 
with a dotted line to guide the eye. The dashed hori- 
zontal line at the top marks the value of the critical r 
for the 3-state Potts model, i.e., the situation at k = 0. 
The dashed curves near the bottom of the plot are the 
results from the perturbative series for small r which we 
briefly discussed in the last section. The Monte Carlo 
data nicely approach these curves for t — 0. 

The curves in Fig. [2] separate the phases with small or- 
der parameter and with (P) /V 1, i.e., the confined and 
the deconfined phases. The phase boundaries depend on 
the inverse mass parameter n, and their behavior is as ex- 
pected for QCD: The intercept of the phase lines with the 
/i-axis shifts to the left with decreasing quark mass (i.e., 
increasing inverse mass parameter k) because a smaller 
chemical potential is sufficient to excite lighter states. 




FIG. 2: Phase diagram as obtained from the maxima of the 
Polyakov loop susceptibility xp- We show results at 4 values 
of the inverse mass parameter k. The dashed curves at the 
bottom are the results of the r expansion and the horizontal 
line marks the critical value of r for the k = case. 

Also the intercept with the r-axis drops with increasing 
K, corresponding to the fact that for quenched QCD (i.e., 
infinite quark mass) the transition temperature is con- 
siderably higher than the crossover temperature of QCD 
with physical quark masses. The mass dependence of the 
phase boundaries thus is as expected for QCD. 

A key problem of the QCD phase diagram is the ques- 
tion about the nature of the various transitions and 
phases. Unless one goes to very high densities where 
more exotic phases exist, two principal phases are ex- 
pected. A phase with conventional matter (confined with 
broken chiral symmetry) and a plasma phase (deconfined 
and chirally symmetric) . In some parameter regions also 
a quarkyonic phase with confinement but restored chi- 
ral symmetry has been discussed. For the transition 
lines a standard scenario is that aX ^ — for physical 
quark masses the finite temperature transition is merely 
a crossover 12] , where different second derivatives of the 
free energy peak at different temperature values [l^] . 
With increasing fi the crossover region narrows and ter- 
minates at a critical endpoint. From the endpoint on 
a first order transition line continues (which at some 
point might hit other transitions to the exotic phases 
mentioned). Alternative scenarios suggest that either no 
critical endpoint appears and crossover type of behavior 
persists also for large fi and low temperature, or that 
even more than one endpoint might exist. For a glimpse 
of the current debate see e.g. [13|. 

While we cannot address questions concerning chiral 
symmetry in the effective center theory, we can ana- 
lyze the type of transitions that take place at the phase 
boundaries. To determine the nature of the transitions 
we used two techniques. We analyzed histograms for the 
distribution of the order parameter and the action near 
the phase boundaries. For a first order transition the 
histograms would display a double peak structure near 
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FIG. 3: Comparison of the phase boundaries obtained from 
the maxima of susceptibility xp a-nd heat capacity C. 

the critical line. In our analysis we found only single 
peaks and thus rule out first order behavior. The second 
approach is a comparison of the normalized susceptibil- 
ities and heat capacities xp/^ a-nd C/V from lattices 
with different volumes. For a first order transition the 
height of the maxima diverges proportional to V, while 
for a continuous transition the divergence is modified by 
a critical exponent. A height which is independent of V 
indicates a smooth crossover. Our analysis shows that 
for the phase boundaries of Figs. [5] and 3 the height is 
independent of V for all volumes we studied. We con- 
clude that the transitions in the effective center theory 
are smooth crossover lines. 

Once the crossover nature is established, one may ask 
how wide the transition region is - similar to the finite 
temperature crossover of QCD at zero density, which is 
20 - 30 MeV wide. In order to get an estimate for the 
width of the crossover region, in Fig. |3] we compare the 
positions of the maxima of the susceptibility xp E^nd of 
the heat capacity C for two values of k. The fact that 
the corresponding lines do not coincide stresses again the 
crossover nature of the transition, and the plot demon- 
strates that the crossover region is rather wide for most 
of the parameter values. Only for small k and fi the lines 
approach each other, anticipating the first order behavior 
which is known for very small k and 0, Q ■ 



Concluding remarks 

In this letter we report on our results for the phase di- 
agram of an effective theory for the center degrees of 
freedom of QCD. Using the fiux representation solves 
the complex phase problem and we develop a general- 
ized worm algorithm for a Monte Carlo calculation in a 
wide range of temperatures and chemical potential. 

The outcome of our analysis is a version of the QCD 
phase diagram when only the center degrees of freedom 
are considered. The phase diagram shares many features 
with the conjectured full QCD phase diagram: The tran- 
sition to the deconfined phase can be driven by both, 
temperature or /i, and the quark mass dependence is 
as expected for QCD. The phase boundaries between a 
phase with only very small center symmetry breaking 
{{P)/V ^ 0) and a phase with {P)/V ~ 1 has a shape 
which is similar to the one conjectured for QCD. For 
all parameter values we studied, the transition is of a 
smooth crossover type and we conclude that center sym- 
metry alone does not provide a mechanism for possible 
first order behavior in the QCD phase diagram. 

Various future research directions may be followed: 
The effective theory can be made more realistic by re- 
placing the Z3 spins by continuous SU(3) valued vari- 
ables (here some work is in progress and also for this 
theory the complex phase problem can be solved by a 
suitable fiux representation). Furthermore it would be 
desirable to take into account also the fermion nature 
of the problem ~ an aspect which is absent in the cur- 
rent effective action. Another interesting direction is of a 
more technical nature: With our effective theory we have 
a reference example of a QCD related system where the 
complex action problem is solved. This reference the- 
ory can and should be used to test the reliability and 
limitations of various techniques for QCD with chemi- 
cal potential, such as reweighting, series expansions or 
complex Langevin methods. 
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